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Abstract 

Magnetorheological fluids (MR) are stable suspensions of magnetizable mi¬ 
croparticles, characterized by the property to change the rheological character¬ 
istics when subjected to the action of magnetic field. Together with another 
class of materials that change their rheological characteristics in the presence of 
an electric field, called electrorheological materials are known in the literature 
as the smart materials or controlled materials. In the absence of a magnetic 
field the particles in MR fluid are dispersed in the base fluid and its flow through 
the apertures is behaves as a Newtonian fluid having a constant shear stress. 
When the magnetic field is applying a MR fluid behavior change, and behaves 
like a Bingham fluid with a variable shear stress. Dynamic response time is 
an important characteristic for determining the performance of MR dampers in 
practical civil engineering applications. The purpose of this paper is to show 
how to use the Optimal Homotopy Asymptotic Method (OHAM) to solve the 
nonlinear differential equation of a modified Bingham model with non-viscous 
exponential damping. Our procedure does not depend upon small parameters 
and provides us with a convenient way to optimally control the convergence 
of the approximate solutions. OHAM is very efficient in practice ensuring a 
very rapid convergence of the solution after only one iteration and with a small 
number of step. 

Keywords: MR dampers, optimal homotopy asymptotic method, auxiliary 
functions, optimal parameters 
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1 Introduction 

Over the past decade, much attention has been given to the MR dampers for 
its attractive characteristics in applications of civil engineering structures including 
earth-quake hazard mitigation, or high strength, insensitivity to contamination, and 
small power requirement. Also MR fluids have attracted considerable interest due 
to their wide range of use in vibration dampers for vehicle suspension systems or 
machinery mounts, their stiffness and damping, characteristics can be adjusted very 
quickly by applying a suitable electric or magnetic field ra-si- The magnetorheo- 
logical response of MR fluids results from the polarization induced in the suspended 
particles by application of an external field. The interaction between the resulting 
induced dipoles causes the particles to form columnar structures, parallel to the 
applied held. These chain-like structures restrict the motion of the fluid, thereby in¬ 
creasing the viscous characteristics of the suspension. The mechanical energy needed 
to yield these chain-like structures increases as the applied held increases resulting in 
a held dependent yield stress [3]. Dynamic constitutive relation of MR fluids is very 
complicated and provided damping force is intrinsically nonlinear, so there is not 
a consistent recognized mechanical model for MR dampers. The mechanical model 
for an MR damper is often established through optimization method according to 
experimental data m-m • At present, there is a variety of dynamic models for MR 
dampers. Some models which are simple cannot effectively simulate nonlinear dy¬ 
namic characteristics of MR dampers. Although there are some models which can 
simulate nonlinear dynamic characteristics, they are established by strong nonlinear 
equations having a lot of parameters which result in complicated numerical calcula¬ 
tion. In civil engineering Bingham model are often used for emulating the dynamic 
behavior of MR dampers. This model is one of the most popular models have been 
widely used with reasonable accuracy and computational cost. The so-called Bing¬ 
ham model includes a variable rigid perfectly plastic element connected in parallel 
to a Newtonian viscosity element. This model assumes that the fluid exhibits shear 
stress proportional to shear rate in the post-yield region and can be expressed as 
0 , 0 , 0 , 0 : 

r = T y (H) ■ sgn( 7 ) + 777 (1) 

where r is the shear stress in the fluid, T y is the yielding shear stress controlled by 
the applied field H, g is the post-yield viscosity independent of the applied magnetic 
field, 7 is the shear strain rate and sgn(-) is the signum function. That is, the fluid 
is in a state of rest and behaves viscoelastically until the shear stress is greater than 
the critical value T y , whereas it moves like a Newtonian fluid when such a critical 
value is exceed. 

Based on this model of the rheological behavior of electrorheological fluids, Stan¬ 
way et al. J5] proposed an idealized mechanical system which consists of a Coulomb 
friction element placed in parallel with a viscous damper. In this model, for nonzero 
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piston velocities x, the force generated by the device is given by the fluid as follows: 


F(t) = f c • sgn(x) + c 0 x + / 0 


(2) 


where f c is the coefficient of the frictional force, which is related to the fluid yield 
stress, Co is the damping coefficient, /o denoting an offset in the force is included to 
account for the nonzero mean observed in the measured force due to the presence 
of the accumulator. We remark that if at any point the velocity of the piston is 
zero, the force generated in the frictional element is equal to the applied force. 
The Bingham model accounts for electro- and magneto- rheological fluid behavior 
beyond the yield point, i.e. for fully developed fluid flow or sufficiently higher shear 
rates. However, it assumes that the fluid remains rigid in the pre-yield region [2]. 
This the Bingham model does not describe the fluid’s elastic properties at small 
deformations and low shear rates which is necessary for dynamic applications. Lee 
and Wereley m and Wang and Gordaninejad Em employed the Herschel-Bulkley 
model to accomodate fluid post-yield shear thinning and shear thickening. In this 
model, the constant post-yield plastic viscosity in the Bingham model is replaced 
with a power law model dependent on the shear strain rate. However due to its 
simplicity, the Bingham model is still very effective, especially in the damper design 
phase [9J. The Herschel-Bulkley model can be expressed by 


r = 


r y (H ) • sgn{ 7 ) + k ( 7 )™ 


sgn( 7 ) 


(3) 


where k is the consistency parameter and rn is fluid behavior index of the magne- 
torheological fluid. For m > 1, Eq. ([3]) represents a shear thinning fluid while shear 
thickening fluids are described by m < 1. For m = 1, the Herschel-Bulkley model re¬ 
duces to the Bingham model m- Zubieta et al. m have proposed field-dependent 
plastic models, for magnetorheological fluid based on the original Bingham plastic 
and Herschel-Bulkley plastic models. In the field dependent Bingham and Herschel- 
Bulkley model, the rheological properties of magnetorheological fluid depend applied 
magnetic held and can be estimated by the following equation [ 12 ] : 

Y = Too + (Y 0 - Too) (2 e~ Ba Ys - e ~ 2Ba Y s ) (4) 

where Y stands for a rheological parameters of magnetorheological fluid such as yield 
stress, post-yield viscosity, fluid consistency and flow behavior index. The value of 
parameter Y tends from the zero applied held value Yq to the saturation value 
Too and a YS is the saturation moment index of the Y parameter, B is the applied 
magnetic density. The value of Yq, Y^ and a YS are determined from experimental 
results using curve fitting method. 


The Bingham body model, presented in Fig. 2 differs from the Bingham model 
(Fig. 1) by the introducing of a spring k. The Bingham body model contains in 
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Figure 1: Bingham model 


parallel three elements that is connecting the elements of St. Venant, Newton and 
the element of Hooke. To a certain value of the applied force f c - static friction force 
of the St. Venant element, only the spring will deform, similarly to the elastic Hooke 
body. If this force is greater than f c the Bingham body will elongate. The rate of 



Figure 2: Bingham body model 

the deformation will be proportional to the difference of the applied force and the 
friction force of the St. Venant element |13j. In theis case the damping force F can 
be expressed as 


p/.x _ / f c sgn(x i) + c 0 ii + f 0 for \F\ > f c . 

\k(x 2 -x 1 ) + fo for \F\ < f c 

where k is the stiffness of the elastic body (Hooke model) and the other parameters 
have the same meaning as in Eq. @. 

An extension of the Bingham model is formulated by Gamota and Filisko [15] . 
This extension describes the electrorheological fluid behavior in the pre-yield and 
post-yield region as well as the yield point. This viscoelastic-plastic model depends 
on connection of the Bingham, Kelvin-Voight body and Hooke body models (Fig. 
3). 
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Figure 3: Extended Bingham model for the Gamota-Filisko model 


The damping force in the Gamota-Filisko model is given by 

{ ki(x 2 - xi) + ci{± 2 - xi) + /o = c 0 x i + f c sgn(x i) + / 0 = 

= k 2 (x 3 - x 2 ) + / 0 for \F\ > f c ( 6 ) 
k(x 2 - xi) + c\_x 2 + /o = k 2 (x 3 - x 2 ) + f 0 for \F\ < f c 

where Co, /o and f c are known from Bingham model ([2]) and the parameters ci, fci 
and ^2 are associated with the fluid’s elastic properties in the pre-yield region. We 
remark that if \F\ < f c then x\ = 0, which means that when the friction force f c 
related with the new stress in the fluid is greater than the damping force F. the 
piston remain motionless. Another view of visco-elastic-plastic properties of MR 
damper behavior is proposed by Li et al. m- In essence the damping force is equal 
to the visco-plastic force, to which besides the friction force connected with the fluid 
shear stress f c , the viscotic force and inertial force contribute, which can be written 
in the form m 

F(t) = f c ■ sgn(x) + cqx + mx (7) 

where Co is a co-factor of viscotic friction and m is the mass of replaced MR fluid 
dependent on the amplitude and frequency of a kinematic excitation applied to the 
piston. 

A discrete element model with similar components, referred to as the BingMax 
model, is reviewed by Makris et al. m It consists of a Maxwell element in parallel 
with a Coulomb friction element as depicted in Fig 4. 

The force F(t) is given by 

F(t) = c - f iie-^-^xir) dr + f c sgn{x{t)) (8) 

Jo 

where g is a parameter meaning the non-viscous damping effect (or the relaxation 
parameter [T 8 ]). The model for the damping force expressed as 

f d (t) = c- f ge-^-^x(r) dr 
Jo 


( 9 ) 
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Figure 4: BingMax model 


was originally proposed by Biot m and later used by several authors in the context 
of dynamics of viscoelastic systems. Eq. ([9|) physically implies that the previous 
time histories of the velocity x contribute to the current damping force and the 
most recent instances of velocity have the highest influence. In the limiting case 


when fi —> oo, the exponential kernel function approaches the Dirac delta function 
5{t). For this special case damping force given by @ reduces to the case of viscous 
damping. For viscoelastic systems, an equation similar to Eq. ® is often associated 
with the stiffness parameter. Cremer and Heckl |20] concluded that: ”of the many 
after-effect functions, that are possible in principle, one - the so-called relaxation 
function - is physically meaningful”. 

Based on the above considerations, in what follows we consider the behavior of 
a magnetorheological damper after Bingham models which in addition contains a 
nonlinear term. 

The motion equation is established in the form: 



+k/3sgn(x(t)) + kf o = 0 


( 10 ) 


where rn. c, fi, k, ak, /3k and kfo are mass, damping relaxation, linear stiffness, 
nonlinear stiffness, coefficient of frictional force and offset of the force, respectively; 
x, x and x are the dynamic responses of the structure (displacement, velocity and 
acceleration). 

The initial conditions are: 



( 11 ) 


The objective of the present paper is to propose an accurate procedure to non¬ 
linear differential equation of the nonlinear Bingham model given by Eq. m, using 
OHAM. A version of the OHAM is applied in this study to derive highly accurate 
analytical expressions of the solutions using only one iteration and a small number 
of steps. Our procedure is independent of the presence of any small or large pa¬ 
rameters, contradistinguishing from other known methods in literature. The main 










MR fluid - OHAM solutions 


7 


advantage of this approach is the control of the convergence of approximate solutions 
in a very rigorous way. A very good agreement was found between our approximate 
solutions and numerical results, which proves that our method is very efficient and 
accurate. 

2 Basic ideas of the optimal homotopy asymptotic method 

Eq. (HOD with initial conditions (I11D can be written in a more general form 

N [s(t)] = 0 (12) 

where IV is a given nonlinear differential operator depending on the unknown func¬ 
tion x(t), subject to the initial conditions 

(13) 

Let xq (t) be an initial approximation of x(t) and L an arbitrary linear operator such 
as 

L[x 0 (t)]= 0, B (x 0 (t), = 0- ( 14 ) 

It should be emphasize that this linear operator L is not unique. 

If p E [0, 1] denotes an embedding parameter and X is an analytic function, 
then we propose to construct a homotopy m-m = 

U [L ( X(t,p )) , H(t, Ci ), N (X(t,p))} , * = 1,2,..., s (15) 

with properties 

n [L (X(t, 0 )) , H(t, a ), N (X(t, 0 ))] = L (X(t, 0 )) = L (x 0 (*)) = 0 (16) 

n [L (X{t, 1 )), Hit, Ci), N (X(t, 1 ))] = Hit, Ci)N (s(t)) = 0 , (17) 

i = 1 , 2,..., s 

where Hit, Ci) ^ 0 , is an arbitrary auxiliary convergence-control function depending 
on variable t and on s arbitrary parameters C\, C 2 , C s unknown now and will 
be determined later. 

Let us consider the function X in the form 

Xit,p) = x 0 (t) +px!(t,Ci). (18) 

By substituting Eq. (fl 8 l) into equation obtained by means of homotopy (fl5l) 

H[L(X(t,p)), H(t,Ci), N(X(t,p))J=0, i = l,...,s (19) 
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and then equating the coefficients of p° and p 1 , we obtain 

U[L(X(t,p)), H(t,Ci), N(X(t,p))] = 

= L(x 0 (t )) + p[L(x 1 (t,C i )) - L(x 0 (t)) - H(t,Ci)N {x 0 (t))] = 0, (20) 

i = 1 , 2 , 

From the last equation, we obtain the governing equation of x’o (t) given by Eq. 
m and the governing equation of x\(t,Ci), i.e.: 

L(x 1 (t,C i )) = H(t,C i )N(x 0 (t)), b(x dXl{t dt Ci) ^j = 0, (21) 

i = 1 , s 

where we find the following expression for the nonlinear operator: 

m 

N {xq( t)) = 22 hi(t)gi(t) (22) 

2—1 

In the Eq. (1221) the functions hj{i) and gi(t), i = 1, m are known and depend 
on the initial approximation xq (t) and also on the nonlinear operator, m. being a 
known integer number. 

In this way, taking into account Eq. (ED, from Eq. (1181) for p = 1, we obtain 
the first-order approximate solution which becomes 


x(t, Ci) = x 0 (t) +xi(t,Ci), i = 1, ..., s (23) 

It should be emphasized that xq( t) and x\ (t, Ci) are governed by the linear Eqs. 
ED and (1211) respectively with initial / boundary conditions that come from the 
original problem. It is known that the general solution of nonhomogeneous linear Eq. 
ED is equal to the sum of general solution of the corresponding homogeneous equa¬ 
tion and of some particular solutions of the nonhomogeneous equation. However, 
the particular solutions are readily selected only in the exceptional cases. 

In what follows we do not solve Eq. ED, but from the theory of differential 
equations, taking into considerations the method of variation of parameters, Cauchy 
method, method of influence function, the operator method and so on, is more con¬ 
venient to consider the unknown function x\(t,Ci), in the form 


xi(t, Cj) = 22 H i(t , hj(t), Cj)gi(t), j = 1, ..., s 


2—1 


tj / r 1 \ dx\{t, Ci . 

B [ X!(t,Ci), - — - I = 0 


( 24 ) 
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where within expression of Hi(t,hj(t),Cj) appear combinations of some functions 
hj, the some terms which are given by the corresponding homogeneous equation and 
the unknown parameters Cj, j = 1, s. In the sum given by Eq. (124|) appear an 
arbitrary number of n of the such terms. 

We have a large freedom to choose the value of n. We cannot demand that 
x\(t,Ci) to be solutions of Eq. (12T1) but x(t,Ci ) given by Eq. (1231) with x\(t,Ci) 
given by Eq. (1241) . are the solutions of the Eq. m ■ This is underlying idea of our 
method. The convergence of the approximate solution x(t, Ci) given by Eq. (1231) 
depends upon the auxiliary functions Hi(t,hi,Cj), j = 1, ..., s. There are many 
possibilities to choose these functions Hi. We try to choose the auxiliary functions 
Hi so that within Eq. ([24]) the term Yli =l Hi(t, hj(t), Cj)gi(t ) be of the same shape 
with the term YliL i hi{t)gi{t) given by Eq. (l22l) . The first-order approximate solu¬ 
tion x(t,Ci) also depend on the parameters Cj, j = 1, ..., s. The values of these 
parameters can be optimally identified via various methods, such as: the least-square 
method, the Galerkin method the collocation method, the Ritz method or minimiz¬ 
ing the square residual error: 


J(C 1 ,C 2 ,...,C S ) = [ R 2 {t,C l ,C 2 ,...,C s ) dt 

Jd 

where the residual R is given by 


(25) 


R(t,C 1 ,C 2 ,...,C s ) = N (x(t,Ci)). 


(26) 


The unknown parameters C\, C 2 , ..., C s can be identified from the conditions: 


dJ _ dJ dJ 

dCh ~ dCh ~ “ dC~s 


(27) 


With these parameters known (called optimal convergence-control parameters), 
the first-order approximate solution given by Eq. (123(1 is well-known. 

It should be emphasized that our procedure contains the auxiliary functions 
Hi(t, fi, Cj), i = 1, ..., m, j = 1, ..., s which provides us with a simple way to adjust 
and control the convergence of the approximate solutions. It is very important to 
properly choose these functions Hi(t, fi, Cj) which appear in the construction in the 
first-order approximation. 


3 Application of OHAM to nonlinear Bingham fluid 
dampers 

In what follows, we apply our procedure to obtain approximate solution of Eqs. 
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cnD and m For this purpose, we introduce the dimensionless variables 



(28) 


such that Eq. m can be expressed as 


x(t) + 2 c • f fie ^ T ^x{r) dr + x(t) + ax 3 (t) + (3sgn(x(t)) + /o = 0 

Jo 


(29) 


where c = g = an d the overdot represents differentiation with respect 

to dimensionless time. The initial conditions (ED become 


x(0) = A, x(0) = vq 


(30) 


Making the transformation 

x(t) = Ae~ xt z(t ) (31) 

where A is an unknown parameter at this moment, Eq. fj2911 can be written as 

z — 2Ai + (A 2 + 2gc + l)z — 2cii 2 ■ f e^ A_M ^ t_r ^ 2 :(r) dr + 

Jo 

+aA 2 e~ 2Xt z 3 + je xt sgn(z - A z) + pp - 2 = 0 (32) 

The initial conditions (1301) become 

~(0) = 1, i(0) = A + ^ (33) 

For the nonlinear differential equation (1321) . we choose the linear operator of the 
form: 


L (. z(t)) = z(t) + u 2 z(t) (34) 

where c o is unknown parameter. 

The nonlinear operator corresponding to Eqs. (l32|) and ([331) is 


N(z(t)) 


—2A z + (A 2 + 2/ic + 1 — uj 2 )z - 
+aA 2 e~ 2Xt z 3 + ^ e xt sgn(z 


2 eg 2 • /' t e (A-A*)(t-T) z ( T ) dT + 

Jo 

- A z) + ^e xt - 2 


(35) 
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The initial approximation zo(t) can be obtained from Eq. (1141) with initial / 
boundary conditions: 


20 (0) = 1, io(0) = - 


co 

T 


(36) 


Eq. (1141) with the linear operator 
(f36l) has the solution 


and with initial / boundary conditions 

(37) 


co 


zo(t) = cos cot —- sinwf 
A 

The nonlinear operator corresponding to Eq. (|22D it holds that 

A(1 + A 2 + ur + 2\ic) — 2/i 2 c(w 2 — A 2 + A^t) 

N (zo(t)) = - w „ . , --cos cot + 


+ 


A (co 2 + (fj, - A) 2 ) 
cj(A 2 + uj 2 — 2 fic — 1)(A 2 + uj 2 + A 2 — 2/rA) + 2^ 2 coo(fo — 2A) 


+ 


A (a ; 2 + (ji- A) 2 ) 

2/r 2 c(a ; 2 - A 2 ) - 2 f icX(co 2 + A 2 - 2\fi) (a _ m)4 
A(w 2 + (/x - A) 2 ) 


sin oot + 


+- 


aA 2 


-2A t 


+ fo e xt + 
A 

3(co 2 4- A 2 ) 3 co(co 2 4- A 2 ) A 2 — 3co 2 

cos cot 4 -—-sin cot 4 -—-cos 3 cot+ 


+ 


A 2 

oo(co 2 — 3A 2 ) 
A3 


A 3 


A 2 


sin Scot 


4/3 


AAl 6 


At 


1 


1 


cos cot -cos 3 cot 4 — cos hoot— 

3 5 


1 1 

— - cos 7 cot 4 — cos 9 cot — 
7 9 


From Eqs. 


4- 


and (1241) . we choose 


z\(t) = C\t 2 cos oot 4- (C^i 2 4- C^t) sin cot + 

co 2 4 - a 2 \ 

(i \/ 2 4~ (uq 4~---7l 4~ 3ACV)f) cos 3utt4~ 


4 -(C 5 t 2 4- C^i) sin 3 cot 4- C 7 ] e 2Xt — C 7 e A * cos wt 


The first-order approximate solution of Eqs. 
obtained from Eqs. (l37l) . (l39l) and (I3T1) : 


and (1301) given by Eq. 


(38) 


(39) 


is 


x{t) = 


+ 


t 2 

co 2 + X 2 ' 


CO 


{C\t 2 + A) cos cot + (C^t" 4 - C^t —-A) sinwf 

A 


e~ xt + 


cos 3cot+ 


Cit 2 + ^ (u 0 + U + x ) A + 3XC 7 j t 

+(C$t 2 4- C§t) sin 3uot + C 7 } e~ 3Xt — C 7 cos cot 


( 40 ) 



















12 


V. Marinca, R.-D. Ene, L. Bereteu 


4 Numerical results 

We illustrate the accuracy of our procedure for the following values of the pa¬ 
rameters: c = 0.1, fi = 20, a = 1, /3 = 0.1, /o = 0.1 A = 5, vo = 0.1. The optimal 
convergence-control parameters are determined by means of the least-square method 
in the three steps as follows: 

For 4 G [0, 7/2] we obtain 


Ci = -8.3835528344, C 2 = 0.4059363155, C 3 = 9.8101433224, 

C A = 12.7300955924, C 5 = -8.3640772984, C 6 = -14.6926248464, 

C 7 = -15.2775231617, A = 0.4221369200, u = 1.1700000000. 

The first-order approximate solution given by Eq. (14011 becomes for this first step: 


x(t) = [(—8.38355283444 2 + 5) cos ut + (0.40593631554 2 + 9.81014332244- 
-13.8580629241) sin cut] e _a4221369200t + [(12.7300955924f 2 - 
-1.02300149694) cos 3 ut + (-8.3640772984f 2 - 14.6926248464) sin 3 ut - 

-15.2775231617] e -1 ' 26641076004 + 15.2775231617 cos cut (41) 

For the second step, when t € [7/2, 7] we obtain 

x(t) = [(—18.71592986704 2 + 1.3931962652) cos ut + (2.58093252414 2 - 
-12.33073302744 - 7.2827523763) sin ut] e -°- 3726878974t + [(16.16871366884 2 - 
-76.62844158384) cos 3 ut + (1.83244632024 2 - 15.48432201434) sin 3 ut - 

-68.7655829044] e - L11 8°636924t + 68.7655829044cos ut (42) 

where u = 1.9481775386, A = 0.3726878974. 

In the last case, when 4 G [7, 10], the first-order approximate solution is 

x(t) = [(0.15215948494 2 + 2.6425574994) cos ut + (-1.42239037514 2 + 
+13.84514477504 - 10.3560770110) sin w4] e -0 ' 1870115667 * + [(-7.99545658504 2 + 
+12.36570702134) cos 3 ut + (1.50736468644 2 - 0.33375012844) sin 3 ut + 

+1.1037184321] g—0-5610347000* _ i.i 0 37184321 cos ut ( 43 ) 

where u = 0.7328908403, A = 0.1870115667. 
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-r(t) 



Figure 5: Comparison between the approximate solution (f^TI) . ll~f2l) and (l43j) and 
numerical solution: -numerical solution; . approximate solution 


In Fig. 5 is plotted a comparison between the first-order approximate solution 
and numerical results. 

It can be seen that the solution obtained by the proposed procedure are nearly 
identical with the numerical solution obtained using a fourth-order Runge-Kutta 
method. 


5 Conclusions 

The Optimal Homotopy Asymptotic Method is employed to propose an analytic 
approximate solutions for the nonlinear Bingham model. Our procedure is valid 
even if the nonlinear differential equation does not contain any small or large pa¬ 
rameters. In construction of the homotopy appear some distinctive concepts as: the 
linear operator, the nonlinear operator, the auxiliary functions Hi(t,Ci) and several 
optimal convergence-control parameters A, to, C±, C* 2 , ... which ensure a fast conver¬ 
gence of the solutions. The example presented in this work, leads to the conclusion 
that the obtained results are of very accurate using only one iteration and three 
steps. The OHAM provides us with a simple and rigorous way to control and adjust 
the convergence of the solution through the auxiliary functions Hi(t,Ci ) involving 
several parameters A, uj, C\, C 2 , ... which are optimally determined. Actually, the 
capital strength of OHAM is its fast convergence, which proves that our procedure 
is very efficient in practice. 
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